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Abstract 

Bioequivalence  is  an  important  area  of  pharmaceutical  research 
containing  many  questions  which  are  not  yet  resolved.  Various  statisti¬ 
cal  approaches  have  been  discussed  in  the  literatures.  We  address  stop¬ 
ping  rules  for  testing  bioequivalence  from  a  decision-theoretic  point  of 
view.  The  numerical  techniques  for  Bayes  sequential  decision  problem 
are  employed  to  obtain  explicit  descriptions  of  the  solutions  of  the 
continuous  time  optimal  stopping  problem  on  bioequivalence. 


Key  words:  Bioequivalence;  Backward  induction;Bayes  risk;  Deci¬ 
sion  theory;  Optimal  stopping;  Sequential  analysis;  Wiener  process. 
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1  Introduction 


Two  or  more  formulations  of  a  drug  are  often  compared  in  a  bioequivalence 
trial.  The  purpose  of  such  a  trial  is  to  determine  whether  alternative  formu¬ 
lations  which  contain  equal  amounts  of  the  same  active  ingredient  give  rise 
to  comparable  blood  levels  or  produce,  in  some  sense,  equivalent  therapeutic 
effects. 

Usually  several  characteristics  of  the  blood  level  concentration-time  curves 
are  considered.  If  a  single  dose  of  the  drug  is  admistered,  then  the  area  un¬ 
der  the  curve  (AUC),  maximum  concentration  (CMAX)  and  the  time  at 
which  the  maximum  concentration  occurs  (TMAX)  all  give  useful  informa¬ 
tion  about  the  the  extent  and  rate  of  absorption  of  the  active  ingredient  of 
the  formulations.  The  distributions  of  these  values  are  usually  not  far  from 
normal.  For  more  discussion  of  the  interpretation  of  relevant  characteristics 
and  design  considerations  associated  with  bioequivalence  trials  see  Metzler 
(1974).  A  lot  of  authors  have  pointed  out  that  a  test  of  the  usual  null  hy¬ 
pothesis  is  inappropriate  since  small  and  clinically  insignificant  differences 
may  be  detected  with  a  large  sample  .  Furthermore,  as  is  always  carefully 
underlined  in  introductory  statistics  courses,  failure  to  reject  the  null  hy¬ 
pothesis  does  not  imply  its  affirmation.  Considerable  controversy  has  arisen 
over  the  appropriateness  of  different  approaches.  For  such  discussions,  the 
reader  is  directed  to  the  articles  by  O’Quigley  and  Baudoin  (1988)  for  general 
approaches  and  Selwyn  et  al.  (1981)  for  the  use  of  the  Bayesian  approach. 

We  address  stopping  rules  for  testing  bioequivalence  from  a  Bayes  se¬ 
quential  decision-theoretic  point  of  view.  Bather  and  Chernoff  (1988)  have 
derived  a  formulation  where  sums  of  successive  observations  of  differences  are 
replaced  by  a  continuous  time  Wiener  process.  Using  several  fundamental 
advantages  of  the  continuous  time  problem  over  the  discrete  time  problem  for 
which  it  is  an  approximation,  they  also  obtained  rough  bounds  and  asymp¬ 
totic  approximations  for  the  solution  of  the  continuous  time  problem.  While 
these  bounds  and  asymptotic  approximations  provide  valuable  insight,  they 
do  not  provide  an  adequate  approximation  to  the  solution. 

In  this  paper  we  will  employ  simple  numerical  techniques  which  are  de¬ 
scribed  in  detail  by  Chernoff  and  Petkau  (1986)  to  obtain  explicit  descriptions 
of  the  solutions  of  this  continuous  time  problem.  The  basic  idea  is  straight¬ 
forward:  the  Wiener  process  is  approximated  by  a  discrete  time  process  and 
backward  induction  is  employed  to  solve  the  optimal  stopping  problem  for 
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this  new  process.  Then  this  solution,  properly  adjusted,  can  be  used  to  ap¬ 
proximate  that  of  the  continuous  time  problem.  This  approximation  can 
then  itself  be  adjusted  further  to  approximate  the  solution  of  the  original 
discrete  time  bioequivalence  problem. 


2  Statement  of  the  Problem 

We  consider  a  trial  with  a  parallel  design  for  comparing  two  formulations,  a 
new  formulation  and  the  standard.  The  design  allows  a  drug  experimenter 
to  terminate  the  program  early  if  the  two  formulations  are  almost  equivalent 
or  far  from  equivalent  and  to  continue  the  trial  otherwise  .  Before  looking  at 
the  sequential  approach  we  recall  the  main  ideas.  Let  /f,  measured  on  some 
scale,  represent  the  true  difference  between  the  two  population  treatment 
means.  In  the  process  we  will  estimate  /x.  As  pointed  out  we  will  be  unable 
to  infer  /x  =  0,  and  even  were  we  able  to  infer  ^  ^  0  this  may  be  of  little 
practical  assistance  if  /x  seems  to  be  close  to  zero. 

The  Bayesian  approach  allows,  in  fact  requires,  explicit  consideration  of 
the  information  available  concerning  the  drug  separate  from  the  current  trial. 
For  example,  the  drug  manager  may  be  quite  sure,  on  the  basis  of  previous 
studies,  that  the  difference  of  the  two  formulations  is  very  small.  On  the  other 
hand,  despite  these  other  studies,  the  bioequivalence  of  these  drugs  may  still 
be  in  doubt,  in  part  perhaps  because  previous  experience  heis  focused  on  a 
rather  narrow  patient  population.  This  information  should  be  used  explicitly 
in  deciding  the  course  of  the  drug’s  clinical  program. 

Consistent  with  the  Bayesian  approach,  the  prior  information  is  quantified 
in  terms  of  a  (prior)  probability  distribution  on  To  be  specific  we  assume 

M  ~  N{fio,  <tI). 

If  Ho  and  (To  are  close  to  zero  then  the  manager’s  prior  assessment  is  that 
the  two  formulations  are  likely  to  be  bioequivalent;  and  large  (Tq  corresponds 
to  a  high  degree  of  uncertainty  regarding  fi-  The  roles  of  /lo  and  (Tq  will  be 
made  clear  in  the  following  development. 

Let  Xi  denote  the  difference  responses  for  z-th  pair  of  patients,  i  = 
1, . . . ,  n.  Assume  the  sequential  random  sample  Xi,....,Xn  are  independent 
N(/x,  xT^),  where  <t^  is  known  (the  unknown  variance  case  is  more  realistic  but 
is  not  conceptually  different). 
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The  posterior  distribution  of  fi  given  Xj, . . . ,  X„  is 
£(/xlXi,...,Xn)  =  X(y„,S„), 


where 

^  _ /^o<7o  ^  +  (Xi  H - +  X„)(T“^ 

(To"  +  n(T-2 

and  the  precision 

=  <To  "  +  n<T"^ 

increases  linearly  in  n. 

So  after  each  observation,  we  need  to  know  n,  the  current  Bayes  estimate 
Yfi  of  fi,  and  its  precision  s„  ;  (ynj^n)  is  the  “state  of  information”  after  the 
n-th  observation.  Then  we  may  decide  to  continue  sampling  or  to  stop.  In 
the  latter  case  we  must  decide  on  whether  or  not  we  have  bioequivalence. 
While  it  will  be  only  approximately  true  in  practice,  we  assume  that  the  cost 
of  the  trial  is  linear  in  the  number  of  pairs  of  patients  in  the  experiment. 
That  is  we  assume  that  the  marginal  sampling  cost  per  pair  is  c.  When  the 
trial  is  stopped,  one  must  decide  to  reject  or  claim  bioequivalence.  The  cost 
of  rejecting  bioequivalence  is  fc,  the  expected  cost  of  having  to  start  over.  In 
the  following  discussion  we  first  consider  the  case  where  the  cost  of  claiming 
bioequivalence  is  /i^. 

Now  let  us  compute  the  posterior  risk  of  stopping  at  stage  n.  We  have 
the  risk  associated  with  stopping  and  deciding  for  or  against  bioequivalence 
plus  the  cost  of  sampling  cn  yielding  s„),  where 

d{y,3)  =  cn  +  min{A:,  =  y,s„  =  s]} 

=  cn  4-  min{fc,  y"  +  s) 

= - 1-  mm{A;,  y"  +  s} - 5-.  (2) 

S  (7q 

The  problem  of  finding  the  Bayes  procedure  for  the  bioequivalence  prob¬ 
lem  has  been  reduced  to  a  standard  stopping  problem  of  the  type  described 
in  Chernoff  (1972). 

Stopping  Problem  :  Let  (Kn,  s„,  n  €  G)  be  a  Gaussian  process  of  inde¬ 
pendent  increments  starting  from  (K„o5^no)»”o  G  G,  with 

C{Yn  -  VmlVm)  =  ^^(0,  Sm  “  5„),  n  >  m,  (3) 
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and  let  the  cost  associated  with  stopping  at  be  given  by  d(V„,s„). 

Find  a  stopping  rule  (a  random  variable  N  taking  on  values  in  G  such  that 
{TV  =  n}  G  B{Yi  :  i  G  G,  no  <  i  <  n})  so  as  to  minimize 

£(d(FA,,SAr)].  (4) 


3  Continuous  time  stopping  problem 

The  continuous  time  stopping  problem  has  a  number  of  fundamental  advan¬ 
tages  over  the  discrete  time  problem  for  which  it  is  an  approximation.  First, 
the  continuous  time  problem  can  be  normalized  so  that  many  of  the  parame¬ 
ters  which  appear  in  the  original  (discrete  time  )  problem  are  eliminated;  thus 
a  single  continuous  time  problem  corresponds  to  an  entire  class  of  discrete 
time  problem.  Second,  the  continuous  time  problem  for  a  Wiener  process 
where  the  cost  associated  with  stopping  depends  only  on  the  stopping  point 
is  related  to  a  problem  in  analysis,  a  free  boundary  problem  involving  the 
heat  equation.  This  relationship  facilitates  obtaining  bounds  and  asymptotic 
approximations  for  the  solution  of  the  continuous  time  problem. 

For  the  bioequivalence  problem,  p  is  regarded  as  a  random  variable,  and 
the  limiting  form  of  the  {Yn,  Sn)  process  is  a  Gaussian  process  of  independent 
increments  Y{s)  in  the  — s  scale  for  5o  >  5  >  5°  ,  where 

E[dY{s)]  =  0,  Var[dY{s)]  =  -ds, 
with  y (so)  =  Ho  at  sq  =  cr^  and 

Note  that  as  time  t  increases  from  0  to  oo,  s  decreases  from  to  0.  Thus 
(—ds)  may  be  thought  of  as  positive.  Hence  a  limiting  form  of  the  bioequiv¬ 
alence  problem  is  a  special  case,  for  G  =  (0,  oo),  of  the  following  continuous 
time  stopping  problem. 

Stopping  problem  :  Given  a  Gaussian  process  {y(s),s  G  G)  of  inde¬ 
pendent  increments  in  the  -s  scale,  with  EdY(s)  =  0,  V'ar[dy(s)]  =  —ds, 
starting  at  y(so)  =  yo,  find  a  stopping  time  S  (5  is  a  random  variable  on  G, 
where  {  5  =  s  }  G  B  {F'(s')  :  so>  s'  >  s})  to  minimize  the  risk 

E[diYiS),S)]. 
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The  continuous  time  version  of  the  bioequivalence  problem  is  essentially 
associated  with  the  cost  function 

ca'i 

d{y,  s)  = - 1-  min{fc,  y*  +  s}, 

s 

after  dropping  the  constant  which  does  not  affect  the  choice  of  the 

optimal  procedure. 

Not  only  is  the  continuous  time  problem  a  limiting  form  of  the  discrete 
time  problem,  but  we  may  regard  the  latter  as  embedded  in  the  continuous 
time  problem  subject  to  the  restriction  that  stopping  may  take  place  only  at 
certain  specified  values  of  s,  i.e.,  s  =  for  integer  values  of  t. 

In  this  continuous  time  framework,  Y  is  regarded  as  a  function  of  s  and  the 
subscript  n  has  been  eliminated  as  an  unnecessary  parameter  which  serves 
only  to  mark  the  possible  stopping  times. 

From  the  point  of  view  of  solving  the  bioequivalence  problem,  certain 
simplifying  transformations  can  be  made.  The  transformation 

y(s*)  =  aF(s) 
s*  =  a^s 

converts  the  F(s)  to  the  F*(s*)  process  which  is  also  a  Gaussian  process  of 
independent  increments  with  F?[dF*(s*)]  =  0,  and 

V’ar[dy“(s*)]  =  a^V’ar[dF(5)]  =  —a^ds  =  —ds*. 

Then  taking  a  =  we  have 

{diy,s)  -  k) 

=  k~^ccr^a^s*~^  +  min{0,  k~^a~^{y*^  +  s*)  —  1} 

=  k~^ccr^s*~^  +  min{0,  y**  +  s*  —  1} 

=  ^4-nun{0,y*’ +  5*  -  1}  (5) 

s* 

Thus  our  problem  may  be  normalized  by  this  transformation  to  that  of  deal¬ 
ing  with  stopping  cost  d*  with  one  sampling  cost  parameter  c*  =  ccr^k~^. 
Now  drop  the  stars,  we  have  a  standard  form  of  optimal  stopping  problem 
with 

=  -  +  min{0,y^ -I- s  -  1}.  (6) 

s 

This  form  involves  only  the  single  parameter  c. 
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4  Numerical  techniques 

Chernoff  and  Petkau  (1986)  have  described  a  number  of  techniques  to  be 
employed  in  obtaining  numerical  descriptions  of  the  solutions  of  the  gen¬ 
eral  optimal  stopping  problem  for  a  zero  drift  Wiener  process  in  the  (y,s) 
scale.  Using  the  same  approach,  we  can  create  a  slightly  modified  program  to 
solve  the  bioequivalence  program.  In  this  section,  we  will  review  the  general 
numerical  techniques  of  Bayes  sequential  decision  problems. 

The  solution  of  a  continuous  time  optimal  stopping  problem  can  be  ex¬ 
pressed  in  terms  of  a  stopping  set  S  and  a  continuation  set  C  =  in  the 
(j/,s)  plane;  that  is,  S  consists  of  stopping  when  (F(s),s)  reaches  as  s 
decreases  from  sq.  This  is  related  to  that  of  a  corresponding  free  boundary 
problem  involving  the  heat  equation.  More  precisely,  that  free  boundary 
problem  is  to  find  («S,6)  so  that 

^byy(y,s)  =  b,{y,s)  for(y,5)eC, 

6(y,s)  =  d{y,s)  for(j/,s)€5, 

by{y,s)  =  dy{y,s)  for(y,s)€95,  (7) 

where  dS  is  the  boundary  of  S.  The  solution  b  of  the  free  boundary  problem 
corresponds  to  the  optimal  risk  d  of  the  stopping  problem,  that  is, 

6(yo,3o)  =  d{yo,so)  =  ^[d(y(S),5)].  (8) 

As  mentioned  before,  the  discrete  time  version  of  the  problem  can  be 
regarded  as  a  special  case  of  the  continuous  time  version  where  stopping 
is  restricted  to  a  limited  subset  of  the  {y,s)  space,  and  hence  the  optimal 
risks  and  related  stopping  sets  are  larger.  In  the  discrete  time  problem,  the 
intervals  between  successive  values  of  s  are  not  equal.  For  convenience  in 
the  numerical  approximation  to  the  solution  of  the  continuous  time  prob¬ 
lem,  we  introduce  another  discrete  time  problem  where  the  successive  val¬ 
ues  of  s  are  equally  spaced.  Moreover,  the  discrete  time  solution  converges 
monotonically  to  the  continuous  time  solution  if  the  set  of  possible  stop¬ 
ping  times  -|-  i  b,  i  =  0, 1, . . .}  increases  and  6—^0.  While  the  value  of 
s  in  the  stopping  times  set  decreases  by  S  between  these  succesive  possible 
stopping  times,  the  process  changes  by  a  normal  deviate  with  mean 
0  and  variance  S;  in  effect,  the  Wiener  process  is  being  approximated  by  a 
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sum  of  independent  normal  random  variables.  At  any  point  {y,s)  where  s 
corresponds  to  a  permissible  stopping  time,  the  choice  between  either  stop¬ 
ping  at  this  point  or  continuing  on  to  next  permissible  stopping  time  and 
proceeding  optimally  thereafter  is  made  on  the  basis  of  which  of  d{y,s)  or 
E[d{Y{s  —  6),s  —  6)|F(s)  =  y]  is  smaller.  Thus,  the  backward  induction 
algorithm  which  yields  the  optimal  solution  to  the  stopping  problem  for  this 
discrete  process  is  specified  by 


d{y,s)  =  d{y,s)  for  s  = 

=  min{d(y,s),£?[d(t/ -f- Z\/?,  s  —  (J)]}  for  s  >  s“,  (9) 

where  Z  represents  a  standard  normal  deviate. 

Note  that  the  evaluation  of  the  expectation  appearing  in  (9)  would  re¬ 
quire  a  numerical  integration  for  which  purpose  the  y-axis  would  have  to  be 
discretized  also.  Thus,  in  practice,  the  backward  induction  is  carried  out  on 
a  grid  of  (y,  s)  points,  each  of  which  is  classified  as  either  a  stopping  or  a 
continuation  point  during  the  course  of  the  computation.  How  would  one  use 
the  results  of  the  backward  induction  algorithm  (9)  to  obtain  approximations 
to  the  boundary  y(s)  of  the  continuation  region  for  the  continuous  time  prob¬ 
lem  ?  Chernoff  (1965)  presents  a  detailed  investigation  of  the  relation  of  the 
discrete  time  normal  problem  to  the  continuous  time  problem.  His  method 
consists  of  simply  adjusting  the  boundary  of  the  optimal  continuation  region 
for  discrete  time  problem;  that  is  the  form 

y5(s)  =  y(s)  ±  0.5826-\/6,  (10) 

where  ys  and  y  represent  the  optimal  boundaries  for  the  discrete  and  contin¬ 
uous  time  versions  and  the  sign  is  determined  so  as  to  make  the  continuation 
region  for  the  continuous  time  version  larger.  This  correction  may  be  used 
to  go  from  the  backward  induction  to  the  continuous  time  version,  and  then 
again  to  go  from  the  latter  to  the  original  discrete  time  problem. 

To  avoid  the  time  consuming  numerical  integration,  the  standard  normal 
deviate  in  (9)  is  replaced  by  a  random  variable  which  is  ±  1,  each  with 
probability  1/2,  leading  to  the  algorithm 

d(y,s)  =  d(y,s)  for  s  =  s®,  (11) 

=  min{d(y,  s),  ^[d{y  -|-  VS,  s  -  6)  +  d{y  -  VS,  s  -  <5)]}for  s  >  s°. 
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We  relate  this  algorithm  to  a  discrete  time  binomial  problem  which  is 
different  from  (9)  which  corresponds  to  a  discrete  time  normal  problem.  For 
the  discrete  binomial  approximation  (11),  the  discretization  of  the  y-axis  is 
necessarily  related  to  the  discretization  of  the  s-axis.  Wherecis  the  Wiener 
process  was  previously  being  approximated  by  the  sum  of  its  increments, 
in  this  simpler  approximation  the  increment  of  the  Wiener  process  is  itself 
replaced  by  a  Bernoulli  random  variable.  While  the  second  moment  of  the 
Bernoulli  variable  is  chosen  to  match  that  of  the  increment  it  is  replacing, 
the  higher  even  moments  do  not  match.  Chernoff  and  Petkau  (1986)  have 
described  another  continuity  correction  for  the  solution  of  the  discrete  time 
version  with  Bernoulli  increments.  Defining 

D{y,s)  =  d{y,3)  -  d(y,s), 

where  d  is  the  optimal  risk  in  the  discrete  time  problem  (  the  function  evalu¬ 
ated  by  the  algorithm(ll))  for  y(s)  6  C  and  close  to  the  y(s),  the  correction 
is  to  use  the  values  Do  and  Dj  of  D  at  yo(s)  and  yi(s),  the  continuation 
points  on  the  grid  which  are  closest  and  second  closest  to  the  stopping  region 
at  the  stopping  time  s,  and  the  continuity  correction  becomes 

J/(^)  =  yo(«)  ±  vy/S, 

where 

A 

^  4Do  -  2Di  * 

The  sign  is  plus  (minus)  when  C  is  above  (below)  S  .  Thus,  by  applying 
two  corrections  we  can  approximate  the  solution  to  the  original  discrete  time 
normal  version  of  our  optimal  stopping  problem.  That  is,  we  calculate  the 
backward  induction  solution  with  the  discrete  time  Bernoulli  process,  use 
(12)  to  approximate  the  solution  to  the  continuous  time  problem  and  end  by 
applying  (10)  to  estimate  the  solution  to  the  original  discrete  time  normal 
version. 


(12) 


5  Implementation 

For  the  bioequivalence  problem,  the  symmetry  of  d{y,s)  about  y  =  0  implies 
that  the  computations  involved  in  the  backward  induction  can  be  confined 
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to  j/  >  0,  that  is,  computing  on  the  grid 

{(j/,s)  :  s  =  s^  +i6,y  =  jVS;i  =  0, 1, . . . , =  0, 1, . . .  (14) 

The  computation  proceeds  in  steps:  At  the  initial  step,  the  values  of  d  are 
assigned  at  all  points  of  the  grid  corresponding  to  the  initial  value  s^,  say 
5^  =  5°  =  0.  We  remark  here  that  in  contrast  to  the  continuous  time 
problem,  the  discrete  time  stopping  problem  under  consideration,  (which  we 
shall  call  the  random  walk  problem),  has  the  property  that  the  continuation 
region  is  truncated;  that  is,  there  exists  an  interval  in  the  s-axis,  [0,  r6]  ,  on 
which  none  of  the  grid  points  will  be  clcissified  as  continuation  points.  In 
fact  for  s  =  6,  d{y,S)  =  d(y,S)  because  d{y,0)  =  cx).  Knowing  d{y,iS)  we 
can  calculate  d(y,  (i  -|- 1)^),  and  it  turns  out  that  for  several  steps,  there  are 
no  continuation  points  and  d{y,iS)  =  d{y,i6)  for  i  =  1,2, ...,r;  Then  for 
z  >  r  +  1 ,  continuation  points  appear. 

In  the  course  of  this  computation  which  yields  the  optimal  risk  for  the 
random  walk  problem,  each  of  the  individual  grid  points  is  classified  as  ei¬ 
ther  a  stopping  point  or  a  continuation  point  for  the  random  walk.  Thus, 
the  continuation  regions  and  their  boundaries  are  determined  and  continuity 
correction  methods  can  be  employed  to  obtain  approximations  to  the  contin¬ 
uous  time  boundaries.  For  accuracy  we  start  with  a  small  step  size  6.  The 
use  of  a  small  grid  spacing  in  a  backward  induction  designed  to  obtain  esti¬ 
mates  for  large  values  of  s  could  require  an  exorbitant  amount  of  computer 
time.  On  the  other  hand,  while  the  use  of  a  large  grid  spacing  may  allow 
the  determination  of  reasonably  good  estimates  at  large  values  of  s,  the  es¬ 
timates  obtained  for  small  values  of  s  would  be  poor.  Thus  the  computation 
is  carried  out  in  stages  or  phases  where  grid  spacings  are  changed  from  one 
phcise  to  the  next. 

The  first  phase  consists  of  starting  at  =  0  and  applying  m,  steps  of 
size  6  for  a  suitably  small  value  of  6.  Then  my,  the  number  of  grid  points 
along  the  y-axis,  must  be  chosen  large  enough  to  contain  all  the  continuation 
points  for  this  first  phase.  In  the  next  phase  we  increase  the  size  of  6  by 
a  factor  of  4  which  automatically  doubles  the  grid  distance  along  the  y- 
axis.  Instead  of  starting  phase  2  at  the  end  of  phase  1  where  s  =  rUgS,  we 
prefer  to  overlap  these  two  phases,  to  give  the  new  coarser  calculation  an 
opportunity  to  adjust,  thereby  avoiding  some  possible  discontinuities  due  to 
the  transition.  Thus  we  have  a  new  6,  four  times  the  original,  and  a  new  s^ 
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between  the  original  and  +  171^6,  and  new  values  of  m,  and  my.  Where 
we  have  overlapping  phases,  we  use  the  finer  grid  to  determine  the  values  of 
the  Bayes  risk  and  optimal  stopping  boundaries  to  be  used  for  publication. 
This  procedure  can  be  repeated  in  successive  phases  of  coarsening  the  grid. 

Referring  to  Figure  1,  we  see  that  there  are  two  boundaries  above  the  y- 
axis.  For  sufficiently  large  values  of  the  constant  c,  the  outer  boundary  turns 
back  toward  the  s-axis.  It  is  possibly  desirable  to  change  S  again  so  that  the 
grid  spacings  become  more  refined  as  the  boundary  gets  close  to  the  s-axis. 
It  is  possible  to  refine  the  grids  by  reducing  <5  by  a  factor  4  when  mo .  mg 
to  the  next  phase.  In  this  case  the  new  will  be  the  last  value  of  s,  i.e., 
-I-  rusS  =  s*.  Now  we  face  a  technical  difficulty.  If  we  label  the  old  and  new 
values  of  (5,  So  and  Sn  =  So/4,  then  the  new  values  of  y  are  i\/^  =  iy/^/l 
and  we  can  not  proceed  because  we  have  not  evaluated  d  at  =  i\/Sll2 
for  the  odd  values  of  i  when  s  =  s*. 

To  overcome  this  difficulty  we  evaluate  d{y,s*)  for  y  =  i\/^/2  for  odd 
values  of  i,  by  replacing  the  last  dichotomous  step  of  ±\/37  by  a  four  valued 
step  with  the  same  mean  0  and  variance  ^o-  In  other  words  if  we  let  y  go  to 


and 


y  ±  with  probability  pi, 


y  ±  with  probability  p2, 


(where  y  =  iy/E^/2  for  odd  i),  then  the  mean  change  E[dY] 
variance  E[dYY  =  So  if 


=  0  and  the 


1 

Pi  +P2  =  - 

Pi  +  9  p2  =  2, 


i.e.  Pi  =  5/16  and  p2  =  3/16.  Thus  for  these  intermediate  values  of  y,  the 
Bayes  risk  at  (y,s*)  will  be  the  minimum  of  d{y,s*)  and 

d{y,  s‘)  =  ^d(y  -  s*  -  So)  +  ^d(y  -|-  -  So) 

+  -  do)  +  ^d(y  +  -  So).  (15) 

Having  calculated  these  values  we  can  now  proceed  with  the  numerical  calcu¬ 
lations  using  the  reduced  value  Sn  of  S.  We  expect  this  technique  would  reveal 
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slight  discontinuities  in  the  estimates  of  both  the  Bayes  risk  and  the  stopping 
boundary  on  moving  from  one  phase  to  the  next.  But  the  experiment  shows 
that  the  jumps  are  so  small  that  we  can  ignore  them. 

6  Numerical  Solutions 

Bather  and  ChernofF  (1988)  have  characterized  the  general  picture  of  the 
solutions  by  studying  the  effect  of  changing  the  standardized  sampling  cost 
parameter  c.  First,  the  optimal  continuation  region  C  will  cover  the  curve 
y  =  ±\/l  —  s  for  0  <  s  <  1.  This  is  because  the  discontinuity  in  first 
derivatives  of  the  stopping  cost  min{0,  +  s  —  1}  implies  a  local  advantage 
in  sampling.  The  advantage  is  of  order  ^|(5s|,  wherecis  the  sampling  cost  is 
of  order  |6s|.  Secondly  C  is  monotone  in  the  sampling  cost  c,  that  is,  Ci  >  C2 
implies  C  <^2.  Third,  every  point  on  the  parabala  y^  +  s  —  1  =  0  belongs 
to  C,  with  the  possible  exception  of  (y,s)  =  (0,1).  In  fact,  there  is  a  definite 
advantage  in  sampling  if  c  <  \J2lte  =  0.484,  t.e.,  (0,1)  6  C  if  c  < 

Fourth,  for  c  >  1,  all  points  (0,s)  lie  in  the  optimal  stopping  set  S  and  for 
0  <  c  <  1,  all  points  (0,s)  with  0  <  s  <  s/c  also  lie  in  S.  Furthermore, 

(y,  s)  €  5  if  c  >  j  and  s  > 

From  the  above  results,  they  have  drawn  the  stopping  boundaries  roughly 
for  c  >  1,  i  <  c  <  y2/7re,  and  sufficiently  small  c. 

We  have  learned  how  the  solutions  would  be  related  to  c,  but  there  is 
no  closed-form  solution  so  far.  While  the  above  results  do  provide  valu¬ 
able  insight,  they  do  not  provide  an  adequate  approximation  to  the  solution. 
Applying  the  previously  discussed  numerical  techniques,  we  explored  the 
sequential  trials  for  a  large  set  of  sampling  cost  parameter  values.  The  nu¬ 
merical  descriptions  of  the  solutions  are  summarized  in  Figure  1  presenting 
the  plots  for  the  continuous  time  version  for  c  =  1.0,  0.5,  0.25,  and  0.1. 

Given  any  fixed  <Tq,  cr^,  fc,  and  c  ,  we  can  calculate  the  optimal  bound¬ 
aries  for  the  original  discrete  problem.  First  we  implement  the  computer 
program  with  the  standardized  sampling  cost  c*  =  c<t*A:“^,  then  apply  (10) 
to  adjust  the  optimal  boundaries  of  the  continuous  time  version  problem. 
As  an  example,  suppose  <Tq  =  5,  <7^  =  20,  A:  =  1,  and  c  =  0.001,  then 


c 

2v/h-  r 
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c*  equals  0.02,  and  the  initial  value  Sq  the  normalized  scale  is  given  by 
•Sq  =  ^~^-So  =  =  5.  The  dotted  curve  in  Figure  2  is  the  optimal  bound¬ 

ary  of  the  continuous  time  problem  in  the  normalized  scale  with  sampling 
cost  0.02.  Two  hundred  optimal  boundary  values  of  the  original  discrete 
problem  are  plotted  and  linked  together  with  solid  line  segments  within  the 
continuation  region  of  the  continuous  time  problem.  Some  of  the  optimal 
outer  and  inner  boundaries  values,  with  notations  and  y'^,  are  listed  in 
Table  1.  The  optimal  sequential  rule  for  this  example  is  to  stop  the  trial  at 
stage  n  if  |y^|  >  or  |y^|  <  and  reject  or  claim  bioequivalence,  and 
continue  the  trial  otherwise. 

One  alternative  model  is  to  consider  that  the  cost  of  claiming  bioequiva¬ 
lence  is  not  but  |^|.  This  leads  to  a  different  stopping  cost. 

E[\y\  I  >^(5)  =  y]  =  5y[Gi(a)  -  ('n(-Qf)],  (16) 


where 

a  =  y  s~^  ,  =  ^{oi)  -t-  a  $(a), 

and  and  $  are  the  density  and  cumulative  distribution  functions  for  the 
standard  normal  distribution. 

Note  that 

Gi(a)  -h  Gi(-a)  =  2  {^(a)  -f  a[$(a)  -  i]}  =  ^i(a). 

We  have 

c^2(y,5)  =  -  +  min  (0  ,  -  1  ),  (17) 

s 

and  we  will  call  the  sequential  optimization  problem  related  to  problem  2, 
and  the  previous  one  problem  1.  It  is  believed  that  the  continuity  regions  of 
this  problem  should  have  shapes  similar  to  those  of  problem  1.  In  particular 
when  s  is  small  both  d(y,  s)  and  d2{y,  s)  are  approximated  by  the  same  term 
c/s,  representing  the  sampling  cost,  and  so  we  expect  similar  behavior  near 
s  =  0. 

The  implementation  of  this  second  version  is  the  same  as  of  the  previous 
one  except  for  replacing  the  cost  function  d(y,  s)  with  ^2(1/,  -s).  Figure  3  shows 
the  shapes  of  the  stopping  boundaries  are  very  similar  for  the  two  versions 
for  four  values  of  c. 
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7  Discussion 


From  the  cost  functions  d{y,s)  and  d2{y,s)  defined  in  (6)  and  (17),  we  see 
the  risk  approaches  infinity  as  s  — >  0.  In  order  to  obtain  better  solutions  for 
s  near  0  ,  Bather  and  Chernoff  (1988)  suggest  the  following  modification  of 
d{y,s)  using  asymptotic  expansion  technology.  Let 

J(a)  =  e-J"'  re^^'dx.  (18) 

Jo 

Note  that  J'(ct)  =  1  —  a  J{oi)  and  both  s~^J{a)  and  $~^J'(a)  are  solutions 
of  the  heat  equation  if  a  =  y  s~^.  Hence,  so  is 

c^-^{J'(a-)  +  ^'(«+)}, 


where 


a_ 


— i —  and  Q!4. 
si 


y  +  1 
1 


Briefly,  the  modification  consists  of  substracting  a  solution  of  the  heat 
equation  from  d  which  does  not  change  the  optimal  solution,  but  makes  the 
derivation  of  asymptotic  expansions  for  the  solution  easier  by  reducing  the 
singularity  due  to  the  term  c/s.  By  subtracting  this  solution  from  d{y,s)  in 
(6)  and  using  J'(q)  =  1  —  aJ{a),  the  new  cost  function  di  is  obtained  : 


di{y,s)  =  min{0,t/^  +  s  -  1}  +  -{a_J(a_)  +  Q+J(a+)  -  1}.  (19) 

s 

Note  the  solutions  of  this  new  cost  function  are  the  same  as  those  of 
d{y,  s).  Thus  we  may  apply  the  numerical  techique  to  di{y,  s)  to  compute  the 
optimal  boundaries  near  (y,s)  =  (1,0).  But  it  is  very  expensive  to  calculate 
the  integral  J(o:)  directly  and  to  apply  the  numerical  approximation  using  di. 
On  the  other  hand  eisymptotic  approximations  to  J{a)  for  small  and  large 
value  of  a  can  be  used  to  provide  asymptotic  expansions  for  the  boundary 
curves  near  the  critical  point  (y,s)  =  (±1,0).  For  small  s;  the  boundary 
behaves  like 


s  2 

Z/x  =  + 

y2  =  1  “  I  -  J  +  "‘25^  (20) 
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where  yi  applies  to  accepting  bioequivalence,  7/2  to  rejection  and  there  are 
symmetric  curves  near  y  =  —1.  The  values  of  ai  and  03  are  |  ~  ^  and 

02  =  I  +  ^  respectively.  Note  that  1  —  |  —  V  ^l^scribes  the  approximate 
behavior  of  the  curve  7/24-5  =  1  near  (7/,  s)  =  (1,0). 

For  problem  2,  Chernoff  and  Bather  also  suggest  that  there  are  almost 
symmetrical  boundary  curves  near  each  critical  point  {y,s)  —  (±1,  0),  at 
y  =  1  ±  0(3^)  and  y  =  —1  ±  O(s^). 

From  the  asymptotic  result  (20),  we  would  expect  to  see  the  two  bound¬ 
aries  close  to  each  other  for  small  values  of  s  for  large  constant  c.  In  order  to 
demonstrate  the  numerical  results  for  small  values  of  s,  we  chose  sufficiently 
small  c  values  and  computed  the  numerical  approximations  for  s  <  1.  Fig¬ 
ures  4  and  5  show  clear  pictures  of  the  behavior  of  the  boundaries  near  the 
critical  points  (y,  s)  =  (±1,0).  We  see,  for  problem  2,  the  estimated  stopping 
boundaries  are  symmetric  around  y  =  1.  The  angles  of  the  curves  get  larger 
and  the  curves  move  forward  as  c  becomes  smaller.  Note  that  a  first  set  of 
numerical  computations  yielded  the  dotted  curves  in  Figure  5  which  did  not 
agree  well  with  asymptotic  expansions  for  s  close  to  zero.  The  continuous 
curves  were  calculated  later,  using  a  smaller  initial  step  size  and  considerably 
more  computer  time.  Even  these  more  refined  calculations  can  stand  some 
improvement  for  s  very  close  to  zero,  where  asympototic  expansions  tend  to 
be  quite  accurate. 

We  are  also  interested  in  how  small  c  must  be  for  the  outer  boundary 
curve  to  never  return  to  y  =  0.  In  general,  we  would  like  to  see  how  the 
inner  and  outer  curves  behave  as  the  sampling  cost  changes. 

We  have  already  learned  from  Bather  and  Chernoff  (1988)  that  the  inner 
curve  and  outer  curve  will  meet  at  (y,  s)  =  (0, 1)  for  c  >  1  for  problem  1.  As 
c  decrejises  to  0,  the  inner  critical  s  value,  the  s  value  where  the  inner  curve 
reaches  the  s-axis,  decrecises  to  0  and  the  outer  critical  s  value  increases  to 
00  .  Table  2  shows  some  of  the  estimated  inner  and  outer  critical  values  for 
problem  1.  Note  that  Bather  and  Chernoff  (1988)  calculate  a  bound  on  the 
inner  critical  s  value  s  >  0.50  for  c  =  0.0554  and  the  numerical  result  is 
5  =  0.5024. 

No  bound  was  calculated  for  that  value  of  c  for  which  the  outer  curve 
never  returns  to  y  =  0.  On  the  other  hand,  the  numerical  calculations 
indicate  that  when  c  =  0.005715,  the  outer  curve  is  still  moving  away  from 
y  =  0  when  s  is  10^*.  We  also  can  get  some  insight  from  Figure  6.  Similarly 
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for  problem  2,  the  results  are  shown  in  Table  3  and  Figure  7.  When  c  >  1 
the  inner  critical  value  and  outer  critical  values  meet  at  s  =  1.57. 
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Table  1 

Estimates  of  stopping  boundary  of  problem  1 
sampling  cost  c  =  0.02 


stage 

n 

outer 

<5  ..out 

■“n  Vn 

inner 

0 

5.0000000 

6.100505 

1 

4.0000000 

5.539054 

2 

3.3333333 

5.085839 

3 

2.8571429 

4.718048 

4 

2.5000000 

4.414596 

5 

2.2222222 

4.167034 

6 

2.0000000 

3.944112 

7 

1.8181818 

3.760117 

8 

1.6666667 

3.596463 

9 

1.5384615 

3.446304 

10 

1.4285713 

3.3169861 

20 

0.8333333 

2.5330751 

30 

0.5882353 

2.1318240 

40 

0.4545455 

1.8871000 

50 

0.3703704 

1.7218550 

60 

0.3125000 

1.6052600 

0.3125000 

0.1656780 

70 

0.2702703 

1.5146050 

0.2702703 

0.2799693 

80 

0.2380952 

1.4439480 

0.2380952 

0.3691222 

90 

0.2127660 

1.3873000 

0.2127660 

0.4412256 

100 

0.1923077 

1.3404170 

0.1923077 

0.5001997 

no 

0.1754386 

1.3018960 

0.1754386 

0.5499748 

120 

0.1612903 

1.2694870 

0.1612903 

0.5918276 

130 

0.1492537 

1.2416210 

0.1492537 

0.6281131 

140 

0.1388889 

1.2176189 

0.1388889 

0.6593656 

150 

0.1298701 

1.1968200 

0.1298701 

0.6866206 

160 

0.1219512 

1.1788880 

0.1219512 

0.7101631 

170 

0.1149425 

1.1629310 

0.1149425 

0.7314749 

180 

0.1086957 

1.1490070 

0.1086957 

0.7504526 

190 

0.1030928 

1.1365730 

0.1030928 

0.7674252 

195 

0.1005025 

1.1308140 

0.1005025 

0.7752527 
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Table  2 

Estimates  of  inner  and  outer  critical  values 
for  bioequivalence  problem  1 


c 

s  -  outer 

wmgm 

1.0000 

1.0000 

bh 

■m 

1.0000 

EIH 

1.0011 

0.50 

0.8973 

1.0164 

0.35 

0.8667 

1.0874 

0.25 

0.7891 

1.2515 

0.15 

0.6778 

2.0060 

0.10 

0.5999 

4.7848 

0.075 

0.5508 

20.059 

0.065 

0.5278 

0.060 

0.5169 

1063.9 

0.058 

0.5102 

19326.5 

0.0575 

0.5083 

172112.6 

0.05725 

0.5082 

4127510.8 

0.0572 

0.5062 

26454535.0 

0.05716 

0.5062 

8290124500.0 

0.05715 

0.5062 

>  2.11671245  (31) 

0.0554 

0.5024 

0.0500 

0.4878 

0.0400 

0.4561 

0.0300 

0.4210 

0.0200 

0.3730 

0.0100 

0.3050 

0.0050 

0.2549 

0.0010 

0.1728 

0.0005 

0.1487 

0.0001 

0.1097 

0.00005 

0.0977 

0.00001 

0.0767 
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Table  3 

Estimates  of  inner  and  outer  critical  values 
for  bioequivalence  problem  2 


c 

s  -  outer 

1.5707 

1.5707 

1.5708 

1.5738 

0.35 

1.6304 

0.25 

1.2841 

1.8043 

0.10 

0.9249 

5.5226 

0.06 

0.7669 

115.62 

0.059 

0.7637 

153.05 

0.056 

0.7491 

512.05 

0.054 

0.7393 

2354.9 

0.053 

0.7339 

10442.4 

0.0525 

0.7319 

41255.5 

0.0523 

0.7301 

100130.5 

0.0521 

0.7299 

413058.6 

0.052 

0.7276 

1401723.5 

0.0519 

0.7280 

18049758.0 

0.0518 

0.7281 

>  2.11671245  (31) 

0.050 

0.7191 

0.010 

0.4150 

0.005 

0.3397 

0.001 

0.2321 

0.0005 

0.2047 

0.0001 

0.1585 
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Figure  2  Boundary  of  original  discrete  problem  1 

C=0.02 


s 
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Figure  4  Stopping  boundary  of  problem  1  --  small  c 
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Figure  5  Stopping  boundary  of  problem  2  --  small  c 
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Figure  6  Inner  and  outer 
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ABSTRACT 


Bioequivalence  is  an  important  area  of  pharmaceutical  research  containing  many  ques¬ 
tions  which  are  not  yet  resolved.  Various  statistical  approaches  have  been  discussed  in  the 
literatures.  We  address  stopping  rules  for  testing  bioequivalence  from  a  decision-theoretic 
point  of  view.  The  numerical  techniques  for  Bayes  sequential  decision  problem  are  em¬ 
ployed  to  obtain  explicit  descriptions  of  the  solutions  of  the  continuous  time  optimal  stop¬ 
ping  problem  on  bioequivalence. 
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